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There are over 100,000 square kilometers of eolian sand dunes and sand 
sheets in the High Plains of the central United States. These land- forms may 
be unstable and may reactivate again as a result of land-use, climate change, or 
natural climatic variability The main goal of this thesis was to develop a model 
that could be used to map an estimate of future dune activity. Multi-temporal 
calibrated Landsats 5 Thematic Mapper (TM) and 7 Enhanced Thematic Map- 
per Plus (ETM-h) NDVI imagery were used in conjunction with the CENTURY 
vegetation model to correlate vegetation cover to climatic variability. This allows 
the creation of a predicted vegetation map which, combined with current wind 
and soil data, was used to create a potential sand transport map for range land 
in the High Plains under drought conditions. 

This study found that Landsat is an excellent tool for mapping current dune 
activity, but the current model of sand transport does not allow for accurate pre- 
dictions of future dune activity. Landsat spatial resolution is necessary to resolve 
small parabolic dunes of eastern Colorado, though MODIS scale imagery could be 
used in conjunction with Landsat data. This study shows that AVHRR scale data 
are far too coarse to use in studying dune activity. Landsat NDVI measurements, 
though imperfect, are readily able to distinguish different levels of vegetation cover 
present on dune crests vs interdune areas. Unfortunately, the available data were 
insufficient to generate a reasonable correlation between NDVI and CENTURY 
output. Further analysis showed that dividing the imagery into larger or smaller 



regions did not significantly impact the quality of the regression, nor does the 
number of Landsat images, within the limited range of images available. This 
study was able to generate a reasonable map of current dune activity for range 
land areas across the High Plains and pointed out some of the problems inherent 
in attempting to predict future activity. 
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Chapter 1 


Introduction 


Extensive eolian sand deposits cover much of the High Plains of the United 
States. These deposits have been mobile at least four times in the past 10,000 years 
(Forman et ah, 1992), and have the potential to remobilize in the future(Ahlbrandt 
et al., 1983; Muhs &: Maat, 1993; Muhs et ah, 1997). Today, these sand dunes 
and sand sheets are largely stabilized by vegetation (Muhs et al., 1996; Muhs k 
Holliday, 2001). Wind speeds are currently high enough that if vegetation were 
to be removed these dunes would likely remobilize. Vegetation in this area is 
highly susceptible to land use practices and periods of drought, such as the dust- 
bowls of the 1930s and 50s. During the 1930s and 50s, droughts, combined with 
poor land management practices, resulted in the reactivation of many of these 
sand deposits(Muhs k Maat, 1993). Because the High Plains are a major source 
of food production for the United States, future eolian activity could have im- 
portant economic and social consequences. This project has attempted to map 
areas of instability in grasslands by combining regional climate data with Land- 
sat 5 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper Plus 
(ETM-b) imagery. 

This thesis is part of a large project studying the High Plains at the Center 
for the Study of Earth from Space (CSES). The long term goal of the project 
has been to map regions of instability and potential reactivation across the High 
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Plains. This larger project has collected three or more images from different years 
over each of thirty-three Landsat scenes in the High Plains and calibrated those 
scenes to apparent reflectance. During the course of the project regions of cloud 
free rangeland, as well as areas of pivot irrigation were hand digitized. 

This project has taken a new approach towards the assessment of dune reac- 
tivation potential. Previous studies have focused on the history of the dune fields 
and regional climate (Forman et al., 2001; Woodhouse k Overpeck, 1998; Muhs 
et al., 1999; Muhs k Maat, 1993). Studies such as Woodhouse and Overpeck(1998) 
have examined the frequency of droughts and dune activity to determine what 
magnitude and duration of drought would be necessary to cause reactivation, as 
well as how likely such a drought is to occur under the current climate regime. 
These studies have been statistically based approaches to estimate the likelyhood 
of dune reactivation in a general sense. This study has taken a physically based 
approach to mapping locations that are likely to reactivate should such a drought 
occur. In the process I have tested the spatial resolution required for such a model, 
as well as the eflPect of temporal resolution. By combining a vegetation growth 
model (CENTURY), climate data, and Landsat Imagery into a single model, this 
project has generated a detailed map of potential dune reactivation for much of 
the High Plains. 

A physically based model of dune reactivation has several advantages and 
disadvantages over previous studies. This model can readily incorporate estima- 
tions of future climate from GCMs. It also generates a map at the same scale as 
Landsat making it potentially useful for regional planning. However, this model 
is based on highly simplified models of the underlying processes and may not ac- 
count for all of the necessary processes. Large-scale statistical studies of the types 
of climatic regimes which have existed in the past inherently account for all phys- 
ical processes. These larger scale studies are only applicable to the future with 
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the assumption that processes will only change as they have historically. Because 
human land use has such a strong effect on the land surface, estimates of future 
dune activity that cannot account for human land use are flawed. Understanding 
the physical processes involved in dune reactivation is important when attempting 
to determine the probable state of the system in a variety of climate and land use 
scenarios. 



Chapter 2 


Background 


2.1 High Plains 

The High Plains extend from South Dakota to Texas, and cover portions 
of Colorado, Kansas, Nebraska, New Mexico, Oklahoma, and Wyoming (fig. 2.1). 
This area receives 200-900mm of rainfall per year, with much of it receiving less 
than 500mm. This classifies much of the area as semi-arid, and makes it highly 
susceptible to both human land use and climatic changes. Over 100,000 square 
kilometers are covered with aeolian sand dunes and sand sheets (fig. 2.1) (Madole, 
1994) which have been active at least four times in the past 10,000 years (Forman 
et al., 1992), with large-scale activity in Nebraska and eastern Colorado occurring 
as recently as 260±60 ^^C yr. BP (Muhs et al., 1997) and evidence suggests 
wide spread activity in the 19th century as well (Muhs & Holliday, 1995). More 
recently, significant regional activity occurred during the dust-bowls of the early 
20th century (Muhs & Maat, 1993). Because a large percentage of the food 
produced in the United States comes from this region, and many people live 
in this area, future changes in land cover are very important both socially and 
economically. In addition, paleo-dunefields are an important record of the regional 
climate over the past lOky, thus a better understanding of the link between 
climate and dune activity may be helpful to paleoclimate studies. 

In the High Plains, vegetation cover is the most important factor controlling 
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High Plains Sandy Soils 



Landsat Image Outlines 

States 

Roads 

Sandy Soils 

Figure 2.1: Sandy soils in the High Plains, based on the STATSGO database and 
Muhs and Holliday (1995). UTM Zone 14 projection, grid is in 100km intervals 
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sand activity. On partially vegetated sand surfaces, the amount of vegetation cover 
has been shown to be the dominant control on sand transport (Buckley, 1987; 
Wolfe & Nickling, 1993; Lancaster & Baas, 1998). Wind speeds in the High Plains 
are often well above the threshold velocity for sand transport, but vegetation cover 
prevents substantial sand movement. In addition the drift potential, a measure 
of the ability of wind in a region to transport sand (Fryberger & Dean, 1979), is 
intermediate to high for the central and northern High Plains (Muhs et ah, 1996). 
In the southern High Plains Muhs and Holiday (2001) concluded that dunes are 
primarily stabilized by vegetation cover. 

Studies have suggested that human induced climate change and natural 
climatic variability are likely to result in future droughts of the magnitude achieved 
in the early 1900s (Woodhouse &: Overpeck, 1998). In addition, human land 
use effects the stability of land cover by influencing vegetation cover and surface 
cohesion. If overgrazing and heavy tilling become common practices, then large 
portions of the sand sheets will have a higher potential to reactivate again in the 
near future. 

Sand dunes in the High Plains are highly susceptible to drought. Vegetation 
in the High Plains is largely water limited (Sala et ah, 1988), thus, droughts are 
periods of high risk for sand activity. Aboveground net primary productivity 
(ANPP) in the High Plains is strongly correlated with precipitation, = 0.9 
(Sala et al, 1988), and ANPP has varied by over ninety percent between wet 
and dry years in portions of the High Plains (Forman et al., 2001). In addition, 
periods of drought favor the expansion of shrubs over grasses (Schlesinger et al., 
1989). Shrubs tend to be more spatially heterogeneous allowing greater erosion to 
occur in the interstices due to lack of cover and a funnelling effect that can locally 
increase the wind speed as air is forced to move around the shrubs. 

The current level of dune activity typically increases from north to south. 
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Figure 2.2: This dune field north east of Kermit, Texas is highly active. 


The most active dunes, such as those outside of Kermit, Texas (fig. 2.2), are 
located in Texas and New Mexico. Eastern Colorado contains large parabolic 
dunes, many of which have moderately active crests and frequent blowouts (fig. 
2.3). The Nebraska sandhills make up the northern extent of this study. This area 
is characterized by large, currently stabilized, transverse dunes with occasional 
blowouts (fig. 2.4). 

2.1.1 Previous Instances of Dune Reactivation 

Sand dunes across the High Plains have been active numerous times in the 
past 10,000 years. Major dune activity is believed to have occurred in the mid 
Holocene, around 8-5 yr BP. Loope et. al (1995) and Stokes and Swinehart 
(1997) present evidence for activity in the Nebraska Sand Hills at this time. For- 
man and Maat (1990) and Forman et. al (1992; 1995) suggest that dune activity 


•SuiSj'ejua moij 

s:^noAvojq [‘bjuij'bu :jiuiq o!^ jjjoAi sjaqou'ej j'eooj uoiSaj aq:; ssojoy -Xij-ejnSaj jnoao 
ni^s s;noMO|q q-eras s^jsaio aunp uo ^nq ‘pazqiq'B:^s qaAiL si 'Baj'e siq:^ jo ^sopq "eqs^jq 
-ajq j'ei:^uao jo qanui jaAoa satinp asaaAsu'ea;^ jo aSBrai ^'Bspu'Bq :p-g amSiq 



•uoiSai siq:^ ui Stiiz-eiS aAisua:^xa 
aq:^ aS'BU'em o:^ pasn si aouaj aqjL 's:^noAio|q Avaj b qqM ‘jaAoo uop^^^aSaA asxeds 
Aq pazqiqB!^s XpSieq -op^jo|oo ii-is^^s'ea m saunp aqoq^jad I'Baid^x -2'Z 9 -mSiq 



8 


9 


in eastern Colorado was also significant during periods of the mid-Holocene. In 
the southern High Plains, Holliday (1989; 1995) found evidence for mid-Holocene 
dune activity in dry valleys. More recently there have been numerous instances of 
eolian activity over much of the High Plains. Studies have shown the past 3,000 
years to be a very active period across the High Plains (Muhs & Holliday, 1995; 
Forman et al., 1992; Loope et al., 1995; Muhs et al., 1996; Muhs et al., 1997; 
Stokes & Swinehart, 1997), the most recent of which probably occurred in the 
19th century (Muhs & Holliday, 1995). 

2.1.2 Effects of Land Use 

Most of the High Plains are extensively impacted by human land use. Range- 
land, used for grazing, and both irrigated and dryland farming are prominent 
through out the region. Some areas have been placed in the Conservation Reserve 
Program (CRP), a program that attempts to increase soil stability by paying 
farmers to plant natural grasses for ten year periods. This study does not look at 
areas which were farmed at any point during the period of study because vegeta- 
tion in these areas responds more directly to variability in human activity then to 
variations in climate. Because most of the remaining grassland is grazed, it was 
impossible to avoid areas which have been impacted by land use patterns, but 
grazed areas are still likely to respond to weather forcings. 

Grazing is a common use for land that is inhospitable to farming. Lands with 
substantial topography, such as sand dunes, or lands that do not have a source 
of water for irrigation and lack the nutrients necessary for productive dryland 
farming are often used primarily for cattle grazing. Variations in grazing can 
affect the vegetation cover substantially, thus weakening the correlation between 
climatic effects and vegetation cover. These effects are noticeable on individual 
Landsat scenes in which fence lines are prominent due to different grazing patterns 




Figure 2.5: The effects of grazing on vegetation cover are visible in this CIR 
image from eastern Colorado. Hard lines between areas with different vegetation 
densities are indicative of recent variations in grazing patterns on opposite sides 
of a fence. 

on opposite sides of the fence (fig. 2.5). 

2'2 Potential Ftiture Climate Change 

Changes in the climate system could substantially affect the amount of mois- 
ture available for vegetation growth as well as wind speeds, and sediment availabil- 
ity in the High Plains (fig. 2.6). Under the current climate, moisture is primarily 
carried to the High Plains by upper level air masses from the Pacific Ocean and 
surface outflow from the Gulf of Mexico (Barry & Chorley, 1987) . Droughts are 
commonly associated with a weaker movement of air from the subtropical Atlantic 
and the Gulf of Mexico (Oglesby, 1991). 

Much effort has been expended attempting to predict future climate. Many 
general circulation models indicate a general warming and drying over the High 
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Figure 2.6: Process model of the effects of climate change on eolian activity in the 
High Plains, after Muhs and Holliday(1995) 
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Plains with a doubling in atmospheric CO 2 (Manabe & Wetherald, 1987; Gregory 
et ah, 1997). Studies have suggested that greenhouse warming may result in a 
climate similar to that of the early to mid-Holocene (Hansen et al., 1988) when 
dune reactivation was more common. In addition to an overall likely decrease 
in moisture availability to the High Plains, an increase in the number of extreme 
events, such as droughts, has been predicted (Gregory et al., 1997; Overpeck et al., 
1990). 

2.3 Sand Transport Theory 

A large amount of effort has also been spent in attempting to understand 
the controls on eolian geomorphology, beginning with Bagnold (1941). Much of 
this work has been spent learning about the physical processes that control sand 
transport. In general, the dominant controls on sand transport are wind speed, 
grain-size, and vegetation cover, though in some locations other factors such as 
soil moisture and cryptobiotic crust can dominate. 

2.3.1 Bare Soil 

The simplest case of sand transport is that which occurs on a flat, bare, 
sandy medium. In this case the primary controls are wind speed and grain- 
size. The seminal work of Bagnold (1941) laid the ground work for the physical 
understanding of sand transport by developing an empirical equation to predict 
the quantity of sand that would be transported in a given wind regime with a given 
grain-size. Further work has gone into understanding the physical processes of lift 
and drag that act on each grain of sand (fig.2.7), and more recently computer 
models have been created that can reasonably simulate the transport of a cloud 
of sand grains (Spies & McEwan, 2000). 
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Figure 2.7. Physical processes controlling sand grain movement, after Pye and 
Tsoar (1990). A sand grain will be lifted off the bed if the difference in pressure 
above and below the grain exceeds the weight of the grain (W). Fi is the lift force, 
Fd is the drag force. Cl and Cs are shape coefficients, p is the density of air, ps is 
the density of the grain, d is the grain diameter, and u is the wind speed. 


Shear Stress vs Wind Speed 

Wind speed is the most important factor in controlling sand transport on 
a bare, sandy surface. Under a laminar flow regime, wind speed varies logorith- 
mically with height above a surface (fig.2.8), thus to understand the force that 
is acting on the surface itself it is necessary to calculate the shear stress (tt*) on 
the surface rather than the wind speed at some height above the surface. This 
is done by accounting for the roughness of the surface and the height at which a 
wind velocity measurement was made with Bagnold’s (1941) wind profile equation 
(2.1). Where z is the measurement height, Zo is the surface roughness length, 
is the wind speed at height 2 , and k is von Karman’s constant, 0.4. Ideally the 
roughness of the surface should be found by fitting Bagnold’s equation through 
wind speed measurements at several heights; if this is not possible a roughness 
length must be estimated. 
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Figure 2.8: Typical wind profile with roughness length Zq. Y-axis is the log of 
height above the surface, after (Wiggs, 1997) 


Threshold Velocity 

Below some wind speed there will not be enough energy to transport sand 
grains. This wind speed is known as the fiuid threshold velocity (u^). The thresh- 
old for a given grain size can be calculated by Bagnold’s (1941) equation (eqn. 
2.2). Where Pa is the density of air, is the density of the solid particle, g is 
the acceleration due to gravity (9.8 p), d is the grain diameter, and A is an em- 
pirical parameter dependent on grain characteristics, roughly 0.1 for sand-sized 
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Figure 2.9. Fluid and impact saltation thresholds as a function of grain size after 
Pye and Tsoar (1990) 


particles. After the initiation of sand transport the effective threshold velocity 
decreases. When saltating sand grains land they impact the grains on the ground. 
At the instant of impact these grains are often lifted a little and thus the amount of 
wind required to maintain saltation is lower than the fluid threshold velocity. The 
velocity needed to maintain saltation is known as the impact threshold velocity 
and is lower than the fluid threshold (fig. 2. 9). 



^ g ^ d 




( 2 , 2 ) 
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Figure 2.10: The effects of vegetation on sand transport, modified from Wolfe and 
Nickling (1993) 

2.3.2 Vegetated Medium 
Effects of Vegetation 

Vegetation exerts four main effects on the potential for sand transport (fig. 
2.10. Vegetation directly covers sand, thereby keeping the wind off of it. Vege- 
tation makes the surface rougher, thus slowing down the wind near the ground. 
Vegetation roots provide a stronger medium than loose sand, and finally, vegeta- 
tion blocks actively moving sand. Many researchers have attempted to quantify 
the effects of vegetation on sand transport. This work can be divided into main 
areas of research, theoretical models which work with generalized roughness ele- 
ments (Gillette & Stockton, 1989; Musick et al., 1996; Raupach et al., 1993) and 
empirical models which have looked at vegetated mediums at specific field sites 
(Wasson k Nanninga, 1986; Stockton k Gillette, 1990; Hagen k Armbrust, 1994; 
Armbrust k Bilbro Jr., 1997; Lancaster k Baas, 1998). A good review of recent 
literature is given in Buckley (1998). 




17 


Theoretical Models 

Many researchers have tried to focus on the more general problem of the 
effects of roughness elements on sand transport. By using simple shapes created 
from wooden blocks, rocks, or other readily mutable objects these researchers 
have hoped to further divide the problem into specific measureable components 
of roughness elements such as height, width, and porosity (Raupach et ah, 1993). 
These techniques are advantageous because if a general solution could be created 
it could then be applied to any type of vegetation, fence, or other barrier. This 
work has had some success with simplified models, but has not scaled well to real 
vegetation due to the complexity of vegetation structure. 

Empirical Models 

The other main approach to predicting sand transport has measured move- 
ment in real world situations or in realistic labratory settings. Empirical ap- 
proaches usually involve correlating wind speed and vegetation cover to sand 
transport in the field (Wasson k Nanninga, 1986; Stockton k Gillette, 1990; 
Hagen k Armbrust, 1994; Armbrust k Bilbro Jr., 1997; Lancaster k Baas, 1998) 
or in a laboratory setting (Buckley, 1987). This approach has the advantage that 
correlations are directly related to real world situations with real plants. It has 
the disadvantage that it is harder to control and results may only be significant 
for the exact sand, vegetation types at the given field site. Laboratory measure- 
ments such as Buckley (1987) are easier to control, but may not have realistic 
vegetation distributions, especially dead vegetation which can play an important 
role in covering sand. Field measurements such as Lancaster and Baas (1998) 
use real-world sites but as a result the researcher has little to no control over the 
vegetation cover or wind speeds. Field work usually yields weaker correlations 
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than laboratory research, but the inability to correlate variables from real field 
sites calls into question the utility of current laboratory measurements which are 
not reproducible in the field. 

2.3.3 Large Scale vs Local Reactivation 

All of these approaches to predicting sand transport focus on local condi- 
tions, but large scale reactivation may be controlled by other factors. Large scale 
reactivation requires that a downwind sink not immediately stop local sand trans- 
port, and local sites can be heavily influenced by upwind vegetation cover and 
sand transport. In addition a significant upwind source may be necessary for long 
term reactivation. Other studies have investigated larger-scale dune reactivation 
but have not developed equations which could reasonably be related to our study. 
Wiggs et al. (1995) correlated their dune activity index to measurements of vege- 
tation cover (R^=0.76) but did not include a measurement of wind speed in their 
study. 

Climate Based Indices 

In an attempt to take a larger scale approach several researchers have created 
climate based regional influences. The most commonly sited of these indices is 
the Dune Mobility Index of Lancaster (1987). The Dune Mobility Index (M) 
relates the percent of time wind speed is above the threshold for sand transport 
on a bare sandy surface {Wt) ^nd the precipitation to potential evapotranspiration 
ratio (P/PE) to dune activity (eqn. 2.3, fig.2.11). 

M = ^ ( 2 . 3 ) 

PE 

Lancaster and Helm (2000) showed that the Dune Mobility Index for the 
previous year is correlated with the amount of sand transport at several field sites 
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in New Mexico and Arizona, but the slope of this relationship varies from 0.0007 
to 0.5 (fig. 2.12). This variation in relationships shows that local factors can sig- 
nificantly alter the relative meaning of the index. While climatic indices give a 
broad picture of the likelyhood of dune activity they do not take account for local- 
ized variations in water resources and vegetation types, making them less useful 
for mapping regions of instability at the scale necessary for land management 
practices. 

2.4 Measuring Vegetation from Landsat 

Landsat data can be used to estimate the percent of vegetation cover at 
each point on the ground. This estimate is far from a perfect measurement. 
Landsat measures the amount of radiation received from within a specific angular 
instantaneous field of view. This includes radiation reflected and emitted by 
the atmosphere as well as radiation reflected and emitted by the ground surface. 
A modest correction was made in this study to account for variations in the 
atmospheric contribution over time, and emitted energy from the earth is, under 
most circumstances, several orders of magnitude less than reflected solar energy 
at wavelengths less than 2.5/im. Thus it is a reasonable assumption that most of 
the variation we see in Landsat data is due to variations in reflected solar energy, 
and in turn related to variability in ground cover. 

The Normalized Difference Vegetation Index (NDVI) (Rouse et al., 1974) is 
a proxy for vegetation cover that can be derived from Landsat data. Vegetation 
typically reflects very strongly in the NIR and weakly in the red portion of the 
spectrum, while soil has a flatter reflectance curve, and senescent vegetation has 
an intermediate curvature (fig. 2. 13). Thus as the amount of vegetation within a 
given pixel increases the total NIR reflectance will increase, and NDVI with it. 
NDVI also removes the effects of variations in topography because it is a ratio 
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Figure 2.11: Dune Mobility Index graph of P/PE vs Wj. Locations with above a 
ratio of 200 are typically highly active, between 100-200 are slightly active, and 
below 100 are generally inactive. After Lancaster (1987) 




Mobility Index for prior calendar year 


Figure 2.12: Correlation between Dune Mobility Index and measured sand trans- 
port at several sites in Arizona and New Mexico, after Lancaster and Helm (2000). 
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Wavelength (micrometers) 


Figure 2.13: AVIRIS reflectance spectra for green vegetation, non-photosynthetic 
vegetation (NPV) and soil. The NDVI for these reflectance curves are approxi- 
mately 0.59, 0.35, and 0.10 for green vegetation, NPV, and soil respectively. (Data 
from (Warner, 2000) 


between two linear band combinations. 

The variation in ground cover may not be limited to variation in the percent 
vegetation cover, thus at sensor radiance will be a function of variables other 
than vegetation cover. NDVI does not account for small variations in soil and 
vegetation reflectance. If soil types vary substantially the contribution of soil 
reflectance to NDVI can vary by as much as 0.3 (Huete et ah, 1985). This means 
that NDVI measurements from different regions may not be directly comparable. 
NDVI is also sensitive to wet soil. Recent rainfall will significantly reduce the 
reflectance of soil and thus alter the NDVI over a wet region. The sensitivity to 
rainfall can make temporal analysis difficult. In addition, NDVI can be sensitive 
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to canopy geometry and non-photosynthetic vegetation (Wessman et ah, 1997). 
Elmore et. al (2000) showed that spectral mixture analysis is a more accurate 
method of determining vegetation cover for temporal analysis than NDVI, but it 
is also requires very careful endmember determination in each scene. 

2.5 Modelling Vegetation with CENTURY 

The CENTURY ecosystem model version 4.0 was used in this study to relate 
vegetation cover to climate parameters. CENTURY was developed at Colorado 
State University to simulate nutrient cycling in temperate grasslands such as the 
High Plains (Parton et ah, 1988). This model has been used heavily in ecosystem 
modeling studies (Gilmanov et al., 1997; Burke et al., 1997; Motavalli et al., 
1994) including climate change scenarios (Parton et al, 1995). Version 4.0 of the 
CENTURY model runs on a monthly time step. The output parameters relevant 
to this study are aboveground live carbon, and standing dead carbon. More 
recently a daily time step version of the CENTURY model (DAYCENT) has been 
developed (Parton et al., 1998), but was not used in this study. 



Chapter 3 


Data and Methods 


3.1 Landsat Images 

This study incorporated over 100 Landsat images from thirty-three path/rows. 
These scenes were collected as part of the larger project at CSES to cover the 
Ogallala aquifer at three separate times. The data were selected in an attempt 
to include a variety of water stress levels. At least one image was collected in a 
drought year, and one in a normal to high precipitation year. Most path/rows 
have one scene from 1985, one from 1988 or 1989 and one from 1996. 1988 and 
1989 were drought years in much of the High Plains. A few path/rows have been 
supplemented with a Landsat 7 scene collected in 1999. This temporal coverage 
allows analysis of the response of vegetation to water stress at each point in the 
High Plains for which we have cloud free coverage. 

3.1.1 Information on Scenes 

The thirty-three path/rows collected for this project cover the area from 
western Texas and southern New Mexico to southern South Dakota and Wyoming 
(fig. 2.1). Generally three to four separate dates for each scene were analyzed, al- 
though eight scenes were available for a detailed examination of eastern Colorado, 
scene 33/32 (tbl. 3.1). 
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Scene 

Image Date 

1 Scene 

Image Date 

Scene 

Image Date 

28/34 

8/1/1988 


30/38 

9/17/1994 


32/31 

5/28/1989 

28/34 

7/22/1996 


30/38 

6/20/1985 


32/32 

8/21/1985 

28/34 

8/25/1985 


30/38 

8/21/1996 


32/32 

7/2/1996 

28/36 

8/1/1988 


31/30 

10/1/1985 


32/32 

9/7/1989 

28/36 

7/22/1996 


31/30 

9/26/1989 


32/33 

9/17/1989 

28/36 

8/25/1985 


31/30 

8/28/1996 


32/33 

8/21/1985 

29/31 

7/15/1985 


31/31 

9/26/1989 


32/33 

7/2/1996 

29/31 

7/26/1989 


31/31 

7/27/1996 


32/34 

9/17/1989 

29/31 

7/29/1996 


31/31 

8/30/1985 


32/34 

8/19/1996 

29/33 

7/15/1985 


31/32 

9/26/1989 


32/34 

9/30/1999 

29/33 

9/25/1988 


31/32 

7/27/1996 


32/34 

7/4/1985 

29/33 

7/29/1996 


31/32 

8/30/1985 


32/38 

9/17/1989 

29/34 

7/29/1996 


31/34 

8/12/1996 


32/38 

8/19/1996 

29/34 

7/31/1985 


31/34 

7/28/1999 


32/38 

7/4/1985 

29/34 

9/9/1988 


31/34 

8/30/1985 


33/30 

8/23/1989 

29/35 

8/14/1996 


31/34 

9/7/1988 


33/30 

8/26/1996 

29/35 

8/24/1988 


31/35 

8/12/1996 


33/30 

8/28/1985 

29/35 

7/31/1985 


31/35 

7/19/1999 


33/31 

8/10/1996 

29/36 

8/14/1996 


31/35 

9/26/1989 


33/31 

8/12/1985 

29/36 

8/24/1988 


31/35 

8/30/1985 


33/31 

9/24/1989 

29/36 

7/31/1985 


31/35 

9/30/1999 


33/32 

7/10/1999 

29/37 

8/14/1996 


31/36 

8/12/1996 


33/32 

8/12/1985 

29/37 

9/25/1988 


31/36 

8/13/1999 


33/32 

8/15/1992 

29/37 

7/31/1985 


31/36 

9/26/1989 


33/32 

6/15/1987 

30/30 

8/23/1985 


31/36 

8/30/1985 


33/32 

9/17/1998 

30/30 

9/3/1989 


31/37 

8/12/1996 


33/32 

9/24/1989 

30/30 

8/5/1996 


31/37 

9/26/1989 


33/32 

8/28/1985 

30/31 

8/15/1988 


31/37 

8/30/1985 


33/32 

10/5/1993 

30/31 

8/23/1985 


31/38 

8/12/1996 


33/33 

8/10/1996 

30/31 

7/4/1996 


31/38 

9/26/1989 


33/33 

9/24/1989 

30/34 

8/15/1988 


31/38 

8/30/1985 


33/33 

8/28/1985 

30/34 

8/5/1996 


32/30 

8/19/1996 


34/30 

8/17/1996 

30/34 

8/7/1985 


32/30 

8/21/1985 


34/30 

8/27/1988 

30/34 

9/7/1999 


32/30 

7/31/1989 


34/30 

7/2/1985 

30/36 

6/20/1985 


32/31 

6/18/1985 


34/31 

7/2/1985 

30/36 

8/21/1996 


32/31 

7/22/1996 


34/31 

9/15/1989 

30/36 

7/31/1994 


32/31 

6/24/1993 


34/31 

8/17/1996 


Table 3.1: Landsat Images used in this study. Landsat scenes are refered to by 
their orbit path/row. 
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Geometric Correction 

As part of the larger project at CSES, all scenes were registered to UTM 
zone 14. Originally all of the 1996 images were warped to vector layers such 
as roadmaps. The next step in the registration process required registering the 
remainder of the images to the 1996 images. To do this, one rough warp, based on 
hand picked ground control points, was performed, followed by a more rigorous 
registration performed with the IMCORR program (Scambos et ah, 1992) which 
typically generates sub-pixel accuracy. Registration was performed with a nearest 
neighbor resampling. 

Radiometric Calibration 

Radiometric calibration is essential for accurate analysis of multi-temporal 
data as well as for meaningful comparisons between scenes. Satellites such as 
Landsat measure radiance at the sensor. Within each spectral band (A), at sensor 
radiance {L\) varies dominantly as a function of surface reflectance (px), solar 
irradiance (Eox), atmospheric transmission down (ta^) and up (tat), the angle 
between the ground surface normal and the incoming solar irradiance {6), and 
path radiance {Lp\) (eqn. 3.1), the subscript A indicates a wavelength dependent 
parameter. Landsat data for the 1996 imagery were converted to radiance using 
the published Landsat gain and offset. These data were then converted to apparent 
reflectance by computing solar irradiance from the sun angle and distance at the 
time of acquisition (from Landsat 7 Data Users Handbook ) , and estimating path 
radiance based on the dark object subtraction method (Chavez, 1996). This 
method does not remove the effect of variations in topography represented by 
cos(6'), these variations are accounted for in the NDVI calculation. 

T _ P*Eox*tx^* ta| cos{6) , ^ 

L>\ - — h LpX 


7T 


( 3 . 1 ) 
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The remainder of the imagery was then cross calibrated to the 1996 re- 
flectance imagery. An advanced pseudo invariant features technique was used to 
calibrate each scene. A spectral angle map (SAM) was created by calculating the 
dot product between the spectral vectors of each pixel from two scenes. From this 
map the pixels with the smallest change in spectral angle were taken to be pseudo 
invariant (their surface reflectance was assumed to be unchanged). For each spec- 
tral band, the pixels from the image to be calibrated were regressed against the 
corresponding pixels from the 1996 imagery to calculate an image to image linear 
calibration. While this calibration process does not entirely remove atmospheric 
effects it does attempt to remove differences in atmospheric effects between dates. 
This method does not account for atmospheric variability between different path 
rows, only between dates for a single path row. 


3.1.2 NDVI 


The Normalized Difference Vegetation Index (NDVI) is a useful measure- 
ment of vegetation that can be derived from Landsat TM and ETM-f- data. NDVI 
is calculated as the difference between Near Infrared (NIR) and red reflectance 
divided by the sum of NIR and red reflectance (eqn.3.2). This method removes 
topographic variability because the cos(0) term from equation 3.1 cancels out 
when you compute a ratio of radiances, this assumes the path radiance has been 
properly removed first. 


NDVI = 


NIR — Red 
NIR + Red 


(3.2) 


3.1.3 Spatial Resolution 

To analyze the importance of spatial resolution, Landsat imagery from east- 
ern Colorado and the Nebraska sand hills were resampled to common lower res- 
olutions. MOD IS and AVHRR are two sensors which could potentially be used 
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in place of Landsat at less cost and with a higher temporal resolution. Landsat 
pixels are 30m across, AVHRR, pixels are 1000m, and MODIS has 250m pixels. 
This means that there are over 1000 Landsat pixels for each AVHRR pixel, and 
just under 70 Landsat pixels for every MODIS pixel. Landsat data were averaged 
together to produce simulated AVHRR and MODIS pixels eastern Colorado and 
Nebraska. 

3.2 Ancillary Data 

3.2.1 Soils Map 

A sandy soils map was created based on the 1:250,000 National Resource 
Conservation Society (NRCS) State Soil Geographic database (STATSGO) using 
the map of eolian sand in the High Plains presented by Muhs and Holliday (1995) 
as a template(fig. 2.1, 3.1). The STATSGO map is divided into a number of 
classes, none of which is clearly eolian sand. Soil types from the STATSGO map 
which appeared to fall within an area designated by Muhs and Holliday (1995) 
as eolian sand were treated as such. The final classes used from the STATSGO 
database are Sandy , Siliceous”, ’’Arenic”, ’’Aridic Paleustalfs” , ’’Ustic Torrip- 
samments , and ’’Ustipsamments” (Mesic and Mixed). Some further hand editing 
was then applied to match individual regions that were on the Muhs and Holliday 
map. This step was necessary to generate a detailed, geographically rectified map 
for model input. 

3.2.2 Rangeland Map 

As part of a larger project at CSES a map of cloud free rangeland across 
the High Plains was hand digitized. This map was created by outlining cloud free 
rangeland in one year for a given path/row then editing that file to exclude any 




Figure 3.1: Muhs and Holliday (1995) map of eolian sand. This map was used 
as a template to generate a map of sandy soils from the STATSGO database (fig. 
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areas that were covered by clouds on other dates. Rangeland was grossly defined 
to cover any area that was not farmed or substantially built-up. This includes 
areas which are grazed as well as other, more natural areas. For simplicity, this 
map includes bed rock, dry lakes, and other features which are clearly non-eolian. 
Ideally the sandy soils map or the rangeland map should eliminate these regions 
because their presence will interfere with the statistics of the correlation between 
NDVI and CENTURY values. Since these odd areas are separated during the clas- 
sification process and cover a relatively small area, they should not substantially 
alter the final results. 

3.2.3 Climate Data 

Precipitation and temperature data were acquired from a large National 
Climate Data Center (NCDC) database (Earthinfo, 2000). This database includes 
the climatic data required to run the CENTURY vegetation model. Only weather 
data that was continuous from 1971 to the final date required for the imagery was 
used. Most data began around 1950 and some records began before 1900. The 
data base is made up of daily measurements, but, because CENTURY is only run 
on a monthly time step, the data were averaged into monthly values. 

3.2.4 Wind Data 

Wind data were taken from the U.S. Air Force Combat Climatology Center’s 
(AFCCC) DATSAV2 Surface Hourly database. These data are gathered at a 
standard 10m height above the ground. Wind measurements are recorded hourly 
in this database. The maximum wind speed in each month was calculated for each 
station. The mean of these maximums was then calculated to generate a simple 
statistic that represents the average maximum wind speed, referred to here as the 
Mean Monthly Maximum. 
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Conversion to Shear Velocity 


To make these measurements useful to the calculation of sand transport 
they were converted to shear velocity. Shear velocity is a measure of the influence 
wind has on the ground surface, and most studies of sand transport relate sand 
transport to wind shear velocity (u*). The roughness around each station was 
unknown, but was assumed to be relatively smooth, a roughness length of 0.0005m 
(zo) was assumed. Roughness length has less effect on wind measurements made 
around 10m, because the term In(^) is insensitive to small perturbations in at 
large values of Most of the effect appears in the lower boundary layer, thus a 
smooth roughness length is a reasonable assumption for these data. Shear velocity 
was then calculated by equation (eqn. 3.3). 


u = 


K*Uz 

In^ 

-2^0 


(3.3) 


Where z is the measurement height (10m), is the wind speed at height z, and 
K is von Karman’s constant, 0.4. 


Spatial Interpolation 

To generate a wind map, these measurements were krigged for each path/row. 
All available wind measurements were searched to find any that existed within 
100km of the path/row boundaries. This t}q)ically yielded nearly 10 stations. 
The Mean Monthly Maximum was computed for each and krigged (Isaaks & Sri- 
vastava, 1989). Kriging was performed with IDL’s Krig2D function using the 
spherical option to fit the semi-variogram, and a large range. The range was 
set at one half the size of the Landsat image. To speed numerical computation, 
kriging was performed at a 300m grid size rather than a 30m grid size. 




0% 1 ocm 


Figure 3.2; Vegetation type maps generated by kriging data points from Paruelo 
and Lauenroth (1996). 


3.2.5 C3 vs C4 Grassland Map 

A map of grass types across the High Plains was generated using published 
values of vegetation types from 73 locations across the Plains (fig. 3.2). Paruelo 
and Lauenroth (1996) compiled a list of vegetation types at 73 locations from 
published papers. They split vegetation types between C3 grass, C4 grass, shrubs, 
forbs, and succulents. All of these categories were krigged across the High Plains 
to generate an approximate map for the entire area. To simplify this further I 
classified succulents and shrubs as C4 and forbs as C3. Though shrubs are often 
C3 and succulents actually use a third photosynthetic pathway (CAM) they were 
lumped in with C4 grasses because they, like C4 grasses, respond more slowly to 
water stress. All five of these categories were normalized to sum to one and the 
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Century Crop Parameter 

% C3 Grass associated with Crop 

CPR 

100 

G4 

75 

G3 

50 

G5 

25 


Table 3.2: Main CENTURY grassland types 


percent of C3 plants was calculated for each grid point (eqn. 3.4). 


%C3 = 


C*3 + forbs 

C3 + forbs + C4 + shrubs + succulents 


(3.4) 


3.3 CENTURY Modeling 

The CENTURY model (Parton et ah, 1988; Metherell, 1992) was used to 
model vegetation abundance across the High Plains. The CENTURY model incor- 
porates vegetation type, soil type, and real weather to model vegetation growth. 
For this study, I used CENTURY Version 4. CENTURY model runs were per- 
formed at all points across the High Plains for which suitable weather data were 
available. 


3.3.1 Parameterization 
C3 vs C4 

The percent C3 and C4 grass was specified in all CENTURY model runs by 
using one of the base parameterizations provided. Built into the CENTURY model 
are 4 classifications of varying C3 vs C4 abundances (tbl. 3.2). For each station 
that was to be simulated in CENTURY the percent C3 was determined from the 
C3 vs C4 grassland map mentioned above and the appropriate parameterization 
was used for that CENTURY run. 
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Weather 

The climate data from the NCDC database were used to parameterize CEN- 
TURY. Daily weather data were converted to monthly averages for all available 
years. In addition, complete monthly means, standard deviations, and skewnesses 
were calculated to allow CENTURY to generate weather stochastically. This 
stochastic weather is used to spin up the model for 2000 years. This allows the 
model to reach an equilibrium with the local climate before trying to match the 
actual weather data. 

Soils 

CENTURY modeling was only performed for sandy regions. All CENTURY 
model runs were parameterized to run with eighty percent sand, ten percent clay, 
and ten percent silt. In reality this parameter varies spatially, but because suffi- 
cient soil maps were not available. Because the only areas I was interested in in 
the end are sandy regions, this assumption is considered reasonable. 

3.3.2 Interpolation 

CENTURY version 4.0 is a point specific model, thus to generate a map 
of modeled vegetation it was necessary to interpolate between point locations. 
CENTURY could only be run at locations for which a viable weather record was 
available, thus there were a limited number of points available. These points 
were interpolated spatially with IDE’s Krig2d function using spherical kriging, 
and a range equal to half of the height of the image. Each path row was krigged 
separately using any CENTURY points located within ten kilometers of the image 
boundary. 
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3.4 CENTURY vs NDVI 

NDVI values were regressed against CENTURY output values in an attempt 
to correlate vegetation cover to weather at each point on the ground. To get a 
better estimate of NDVI, large areas of similar vegetation cover were averaged 
together, these areas were determined by a classification process described below. 
NDVI is influenced by the amount of non-photosynthetic vegetation as well as 
green vegetation. To account for the influence of non-photosynthetic vegetation 
half of the CENTURY output variable stded (standing dead carbon) was added 
to the CENTURY value for aglivc (above ground live carbon). This final value 
was regressed against NDVI. 

3.4.1 Temporal Classification 

To find areas of similar vegetation, NDVI time sequences were classified 
with the ISO-data classifier. Details of the ISO-data classifier can be found in 
Tou and Gonzalez (1974). All available NDVI images for a given path row were 
concatenated together and masked to eliminate areas that contain clouds, non- 
sandy soils, or have been heavily influence by humans. The resulting image was 
classified on the basis of the time variance of NDVI at each pixel. Each class 
in the resulting image should represent a group of pixels in which the vegetation 
responds similarly to recent weather. To investigate the importance of the number 
of classes used, the classifier was initialized with 5, 10, 15, 20, 25, 30, 35, and 40 
classes in separate runs. Each class was required to have at least 1000 pixels in it. 
The classifier was allowed to join and separate classes to find the optimal clusters, 
and the classifier was run for three iterations. 



35 


3.4.2 Correlation 

From the maps generated by interpolating CENTURY values for agliv and 
stded, it was possible to correlate NDVl with CENTURY both spatially and tem- 
porally. This regression model tries to find the correlation for each class of ground 
cover mapped by the ISO-data classifier. Every reference to correlation is on a 
class by class basis. To look for the combination of spatial and temporal relation- 
ships, NDVI values were regressed against every corresponding CENTURY value 
within a class. This form of correlation has a data point for every pixel in every 
available image. To look for temporal correlation, all of the NDVI values in each 
ISO-data class were averaged together, as were all of the CENTURY values; this 
generates a single point for each date of available imagery. Each of these averages 
was computed for each date of Landsat imagery available and the results were 
regressed with ordinary least squares regression. Because most scenes only have 
three years of imagery the temporal regression is severely limited in the number of 
data points, therefore, the combined spatial and temporal correlation was used for 
the model output. The temporal correlation was used to investigate the impact 
of the number of classes and the number of images used in the regression. 

3.5 Future Vegetation 

A future vegetation map was generated based on the correllation between 
NDVI and CENTURY for each class. A potential extreme drought CENTURY 
value (30g Carbon) was used with the regression developed previously to estimate 
NDVI for each class during an extreme drought. This value was then converted 
to vegetation cover by assuming that the minimum and maximum two percent 
of NDVI values from the original Landsat imagery correspond to 0% and 100% 
vegetation cover respectively and interpolating between those values. Because 
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the correlation in some of the classes was negative or had an extremely low 
value a statistical method was used as a fall back estimate of vegetation cover. 
If the correlation was deemed spurious, the average NDVI value for a class was 
calculated and divided in half to represent drought conditions. 

3.6 Sand Transport 

A potential sand transport map was generated based on the future vegeta- 
tion map and the wind map. At each pixel the equation derived by Lancaster and 
Baas(1998) (eqn. 3.5) was used to estimate potential sand transport. 

q = 300 * (u* — u*t)^ * g-25*A ^2 

where : 

q = Total sand flux 

u* = Wind shear velocity 
«*t = Threshold shear velocity 

A = Vegetation cover 

3.7 Model Overview 

The entire sand transport model is illustrated in figure 3.3. Precipitation 
and temperature data are used in conjunction with maps of vegetation type and 
soil type as inputs to the CENTURY vegetation model to predict potential veg- 
etation across the High Plains. Separately, Landsat measurements of vegetation 
(NDVI) are classifled to map areas of similar ground type. A time series of NDVI 
measurements for each of these ground types is correlated with a time series of 
potential vegetation predicted by CENTURY. These correlations are used to pre- 
dict vegetation cover for each class of ground type in a drought scenario. This 



37 



Figure 3.3: Sand transport model overview. Basic model inputs are on the left, 
the process is described in detail in the text. 


map of future vegetation is then combined with a map of present wind speeds to 
predict potential sand transport across the High Plains. 








Chapter 4 


Results 


4.1 Areas of Apparent Instability 

According to this model, the largest areas of potential instability are the 
Nebraska Sand Hills, and portions of eastern Colorado, New Mexico and Texas 
(fig. 4.1). Isolated locations through out the High Plains also appear to be at ri s k, 
but no place else has as large an area which appears likely to reactivate. 

4.1.1 Arecis with Clear Dune Landforms 

Many of the areas that this model mapped as having a high potential for 
reactivation exhibit clear eolian land forms. One would expect that locations that 
have experienced recent sand activity would be the most likely to reactivate in 
the near future. The region south of the South Platte river in Eastern Colorado 
has many clear parabolic dune forms (fig 4.2), evidence that this region has been 
active many times in the past. Areas here that are mapped with an extremely high 
reactivation potential are commonly semi-active today with numerous blowout 
features. For the areas in Eastern Colorado multiple model runs were analyzed, 
in all cases these areas were mapped with a high reactivation potential. 

Other areas that were mapped with a high potential for reactivation include 
the Nebraska Sand Hills, and portions of western Texas and eastern New Mexico. 
The Nebraska Sand Hills are large transverse dunes (fig. 2.4), and many of the 









Sand Dune 
Activation Potential 



Figure 4.1: Potential sand dune activity map for the entire High Plains. Color 
scale indicates potential for sand dune activity. White designates areas with 
negligible activation potential and regions that were not mapped by this study. 
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Figure 4.2: Parabolic dunes are common in Eastern Colorado. This area, South 
East of the South Platte river, shows a various levels of dune activity. The 
parabolic shapes in this images are only moderately elongated. 


locations in Texas are currently active sand dunes (fig. 2.2). 

4.1.2 Areas with Clear Non-Eolian Landforms 

Some areas that were mapped as having a high potential for reactivation 
actually exhibit non-eolian land forms. One area in particular in Eastern Colorado 
shows clear fluvial erosional features (fig. 4.3). This area fell with in scene 33/32, 
thus eight images were available for the model. The sensitivity analysis used to 
study the effect that the number of available Landsat images had on the model 
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(a) (b) 

Figure 4.3. This section of Northern Colorado has sections mapped with a high 
reactivation potential (red and yellow, [a]). Examination of a CIR Image of the 
areas shows clear fluvial features [b]. 


showed mixed results for this area, but the bright stream bed was always mapped 
with a high potential for reactivation. 

4.2 NDVI vs CENTURY Correlation 

In general the correlations between NDVI and CENTURY modeled vegeta- 
tion cover were poor. The coefl&cient of determination was calculated for each class 
in each scene, for all points, and for temporal correlation. The number of classes 
used by the model was allowed to vary from 5-50, and the coefficient of determi- 
nation was calculated to analyze the effects of this modification in the number of 
classes. For scene 33/32 the coefficient of determination was also calculated with 
varying numbers and dates of input scenes. 
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Effect of the Number of Classes on Correlation 
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Figure 4.4: Coefficient of determination from temporal regression of NDVI to 
CENTURY plotted against the number of classes used to seperate different ground 
types. 


4.2.1 Sensitivity to the Number of Classes 

There was no apparent trend in the degree of correlation between NDVI 
values and CENTURY values due to the number of classes in the model. The 
coefficient of determination was graphed against the number of classes for all 
points, and for the temporal correlation alone (fig. 4.4). 

4.2.2 Sensitivity to Temporal Resolution 

In scene 33/32 it was possible to analyze the effect that the number of 
available Landsat images has on the degree of correlation between NDVI and 
CENTURY. Random images were chosen for 3-7 image subsets. There was an 
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Figure 4.5: Coefficient of determination from regression of NDVI vs CENTURY, 
plotted against the number of images used in the model. 


apparent decreasing relationship between the amount of variation explained by 
CENTURY with an increasing number of input images (fig. 4.5). This is probably 
due to the fact that with only three images it is very easy to get a spurious 
correlation. Examination of the p- values from the F-test shows no correlation 
between quality of fit and the number of images (fig. 4.6). 

4.3 Spatial Resolution 

Landsat images which were spatially averaged to simulate AVHRR and 
MODIS imagery show a substantial loss in detail. In eastern Colorado dune 
crests are blurred in MODIS imagery, but crop circles are still resolved(fig. 4.7b). 
AVHRR pixels sizes do not definitively distinguish between sparse agriculture 
and open rangeland, and dunes are completely obscured (fig. 4.7c). In Nebraska 
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Effect of the Number of Images on Correlation 



Niunber of Images 

Figure 4.6: The average F-statistic p-value for the temporal correlation of NDVI 
and CENTURY, plotted against the number of images used in the model. 


MODIS data are able to clearly distinguish major transverse dune features and 
separate the large, heavily vegetated interdune areas from the dune crests (fig. 
4.8b). In this same area AVHRR is unable to resolve the dune crests, thus large, 
potentially unstable dunes are averaged together with heavily vegetated interdune 
areas (fig. 4.8c). Note that although it is possible to resolve the major transverse 
dunes in Nebraska with MODIS, it is not always possible to differentiate lakes 
from dune crests. However, this would be easy to do in a CIR image. It is impor- 
tant to note that these represent best case scenarios for both MODIS and AVHRR 
with perfect sampling, and extremely low noise because they are the result of av- 
eraging Landsat pixels together, real AVHRR and MODIS data would have much 
more noise, and delineation of features that are only one pixel wide would become 
harder. 
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Figure 4.7: Low spatial resolution imagery in eastern Colorado shows the limita- 
tions of the MODIS and AVHRR sensors, a) is the original Landsat image, b) is 
the Landsat image resampled to MODIS (250m) pixel size, and c) is the Landsat 
image resampled to AVHRR (1km) pixel size. 



Figure 4.8: Low spatial resolution imagery in the Nebraska sand hills shows the 
limitations of the MODIS and AVHRR sensors, a) is the original Landsat image, 
b) is the Landsat image resampled to MODIS (250m) pixel size, and c) is the 
Landsat image resampled to AVHRR (1km) pixel size. 




Chapter 5 


Discussion 


5.1 Vegetation Measurement 

Neither NDVI nor CENTURY represent a perfect measurement of vege- 
tation. NDVI values within rangeland were assumed to be unaffected by any 
processes other than vegetation growth due to weather and localized conditions 
such as topography, which are assumed to be time invariant. This is clearly not 
the case as much of the area is grazed in a time variant fashion, and recent rain can 
strongly effect NDVI values by darkening the soil. CENTURY values are modeled 
vegetation estimates which are completely dependent on the parameterization of 
the model. 

5.1.1 Grazing 

Areas that are heavily grazed pose a problem to this study. From Landsat 
data, it is impossible to readily distinguish areas which have been grazed from 
areas which have not. Heavy grazing during a wet year and a lack of grazing 
during a dry year may cause the vegetation at a given locality to apparently 
respond inversely to what would be expected under climatic forcing alone. The 
effects of grazing are clearly present in the images. Occasionally fence lines are 
visible due to differential grazing patterns (fig. 2.5). A tighter integration of land 
use studies could possibly alleviate some of these problems as it is possible to 
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include grazing in the CENTURY model. 

5.1.2 NDVI 

NDVI is not necessarily directly correlated with vegetation cover. NDVI 
is based on the variation in apparent reflectance between the red portion of the 
spectrum and the NIR portion (flg. 2.13). Because vegetation has a much higher 
difference in reflectance between the NIR and red than soil, a simple measure of 
this difference is often highly correlated with the amount of vegetation cover in a 
pixel; this is especially true for arid environments. Numerous studies have shown 
that Leaf Area Index (LAI) varies linearly with NDVI for LAI values less then 2-3, 
depending on leaf angle, and then approaches some maximum value asymptoti- 
cally (Tucker, 1979; Asrar et ah, 1984; Neman! & Running, 1989; Carlson et ah, 
1990; Carlson & Ripley, 1997). And LAI and fractional vegetation cover are also 
strongly related, as illustrated by Carlson and Ripley (1997) . However, different 
soil types can have different NDVI values, and within a specific soil type NDVI 
can vary substantially due to moisture content. Because this model is only de- 
signed for sandy soils the variability of NDVI due to varying soil type is probably 
not very large, but it could still be significant. The variability in NDVI due to 
soil moisture is much harder to account for and may be a large source of error 
in this study. Soil moisture decreases NIR reflectance, thus lowering NDVI for a 
region. Summer thunderstorms are common in the High Plains, and can be highly 
localized and ephemeral, making it difficult to account for, especially at the scale 
of Landsat. Given a map of rainfall prior to each Landsat scene a mask could be 
generated to throw out any NDVI values which are likely to have been effected. 



5.1.3 


CENTURY 


CENTURY model output may not be representative of vegetation cover 
at the time of a Landsat overflight. The CENTURY model has been show to 
be accurate within 25% in areas of semi-arid grasslands (Parton et al., 1993), 
but it may be prone to substantial errors due to imprecise parameter estimates. 
Because this study has attempted to extend CENTURY measurements over a 
broad area certain weaknesses are inherent to our use of it. It was not possible 
to verify every measurement location to determine exact soil type, nutrient con- 
ditions, or vegetation types, thus our estimations of these quantities stem from 
a general sandy soil estimate, and a vegetation map generated from point values 
spread across the High Plains. In addition, CENTURY output is on a monthly 
basis. Landsat measures vegetation at a single point in time. Recent rainfall may 
have allowed substantial vegetation growth in the days leading up to a Landsat 
overflight. Herbaceous vegetation can respond very rapidly to rainfall events and 
change on this time scale is not predicted with CENTURY’S monthly time step. 

Vegetation type 

Vegetation type is a very important parameter in CENTURY which may 
be improperly estimated in this study. Vegetation type determines how efficiently 
plants use water and nutrients. C4 plants tend to be much better at conserving 
water than C3 plants due to their specialized photosynthesis process. To attempt 
to include this variability in the CENTURY model, a map of C3 vs C4 vegetation 
type was created from points compiled by Paruelo and Lauenroth(1996) (fig. 3.2). 
The points compiled in this paper actually refer to C3 grasses, C4 grasses, forbs, 
shrubs, and succulents. CENTURY does not explicitly model forbs or succulents, 
though it can model shrubs. Because of the added complexity involved in modeling 
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shrubs, this model further simplifies everything into C3 or C4 grasses. Forbs 
typically use the C3 photosynthesis pathway, thus they were lumped into C3. 
Shrubs commonly use the C3 pathway, but respond more slowly to changes in 
precipitation so they were lumped as C4. Succulents use a third photosynthetic 
mechanism classified as CAM. Because the CAM pathway is more similar to the 
C4 pathway, succulents were also lumped into the C4 group. This lumping is 
partially justified by the fact that the based on the maps in figure 3.2, forbs 
and succulents appear to play a relatively minor role, and even shrubs are only 
dominant in a few locations. 

This classification into C3 and C4 vegetation type throws out some of the 
information available; a better model may result from more detailed vegetation 
classification. Forbs typically respond more rapidly to rainfall, and have a strong 
infiuence on NDVI, so their inclusion as C3 may underestimate changes in NDVI 
over time, especially if there has been very recent rainfall. Succulents and shrubs, 
on the other hand, typically respond much more slowly to variations in climate 
than do C4 grasses, thus lumping them in with the C4 class may increase CEN- 
TURY’S response time for vegetation due to weather forcings. In addition, the 
spatial estimate of vegetation type is based on the assumption that vegetation 
type varies smoothly across the High Plains. This is certainly not the case as 
vegetation type often varies substantially due to localized variability in water and 
nutrient availability, not to mention human infiuence. 

Soil Type 

CENTURY also provides for detailed soil parameterization. Because we 
were limiting this model to areas of sandy soil, all CENTURY runs were performed 
with 80% sand, 10% clay, and 10% silt. The STATSGO Database used to derive 
areas of sandy soil includes more detailed information on soil types, but there is 
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Figure 5.1: Comparison of a high spatial resolution image (a) (Ikonos 4m CIR 
composite) to a classification developed firom a time series of Landsat (30m) NDVI 
measurements (b). The yellow class indicates relatively stable ground, the red 
class represents areas that are partially stabilized, and the gold class represents 
areas that are nearly active today. 


no clear mapping between their soil types and CENTURY’S soil types. This is an 
area that could potentially be improved in the future with a careful interpretation 
of the STATSGO database. 

5.2 Classification 

The classification algorithm appears to produce a useful division of the land- 
scape. The resulting classes visually appear to be consistent in their ground type. 
Comparing a close up of the area around Roggen, Colorado with a Landsat scale 
classification image shows that specific classes are consistently related to specific 
ground types (fig. 5.1). 

In order to examine the effect of the number of classes on the relationship 
between CENTURY and NDVI, a sensitivity study was performed. The classifier 
was parameterized to create between 5 and 50 classes in increments of 5. It was 
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still permitted to join classes together, but it would only do so if classes were 
extremely similar. The relationship between CENTURY” and NDVI was then 
recomputed for each classification scheme. 

Coefficients of correlation (R^) were compared with the number of classes 
used to generate them. Figure 4.4 shows a plot of R^, normalized by the number 
of points used in each regression, against the number of classes. Based on the 
T-Means test here is no significant difference between the R^ values computed 
with greater then or equal to 20 classes and R^ values computed with less than 
20 classes. This indicates that the number of classes used when comparing CEN- 
TURY to NDVI is not very important. From analysis of the maps generated 
by this project, the number of classes used for each Landsat scene (from 5 to 15 
classes, depending on the scene statistics) appears to be appropriate because dune 
crests and large sandy areas are clearly differentiated from interdune and other 
heavily vegetated areas. 

5.3 Landsat Temporal Resolution 

One of the greatest weaknesses of the present model could potentially be 
alleviated with higher temporal resolution satellite coverage. If the poor corre- 
lation between Landsat measurements of NDVI and CENTURY model output is 
purely due to random noise in the samples, then a greater number of samples 
may be the only way to improve the statistical correlation. For most of the scenes 
analyzed in this study only 3 years of imagery were available to correlate NDVI 
with CENTURY. In a temporal correlation, using three points to define a line 
leads to statistically poor correlations. For scene 33/32 (Northeastern Colorado) 

8 dates were available ana sufficiently cloud free to be analyzed. However, even 
8 separate images did not generate a good correlation (Fig. 4.5). This may be 
due to other weaknesses in the model, or it may be that a very high temporal 
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resolution data set, such as is available from MODIS, would be necessary. It does 
indicate that doubling the number of Landsat images used from 3 to 6 is unlikely 
to substantially improve the results. 

5.4 Landsat Spatial Resolution 

Landsat spatial resolution appears necessary to accurately map sand dunes 
in the High Plains. The two other sensors that would allow higher temporal 
resolution NDVI imagery are MODIS and AVHRR. MODIS has a 250m pixel 
size, and AVHRR has a 1000m pixel size. Simple resampling results (figs 4.7 
and 4.8) suggest that AVHRR data are largely useless for mapping sand dune 
reactivation, and that MODIS data are potentially very useful. 

MODIS data can resolve most important features. It can not resolve in- 
dividual dune crests in eastern Colorado, but dune crests tend to occur in close 
proximity, and the combination of Landsat data with MODIS data may prove ad- 
equate to study this area. Dunes in the Nebraska Sand Hills are readily resolved 
by MODIS, thus MODIS data is considered adequate to study these large linear 
dunes. 

It is impossible to resolve important features such as the major linear dune 
crests in the Nebraska sand hills with AVHRR data. AVHRR would be an enor- 
mously useful resource because it has an extremely high temporal resolution and 
a long history. Because the dunes can not be resolved (fig. 4.8) it is impossible to 
determine how much signal is coming from dune crests and how much is coming 
from healthy vegetation in interdune aresis. As a result there is no way of analyz- 
ing the vegetation on dune crests exclusively, even if Landsat data are used ahead 
of time to make an accurate map of the area. 
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5.5 Wind Speed 

5.5.1 Single Index vs Temporal Integration 

Wind measurements used in this study are meant to be an index of the 
general windfield across the High Plains, but wind varies substantially through 
time. Each point measurement used in this study contains a reasonably complete 
temporal sequence of hourly measurements. Ideally all of this data should be used 
and sand transport should be integrated under the time sequence of wind speed. 
Because it would also be necessary to generate a map of wind speed at each time 
step this process is not practical, so some approximation must be used. 

Sand is dominantly transported by only the strongest winds, thus it makes 
sense to look at the average high wind speed for each location. I computed the 
mean monthly maximum wind speed as such an approximation. Other approxi- 
mations might include looking at the top 2% cutoff, or incorporating the percent 
of time the wind stays above some threshold value. The top one and two percent 
cutoff values were analyzed and determined to be highly correlated with the mean 
monthly maximum (fig. 5.2), though the mean monthly maximum was always 
higher than either of these values. Because they are so well correlated it was 
determined not to be important which of these approximations was used and the 
original mean monthly maximum value was used. 

Wind often changes seasonally, but the present single index should still be 
reasonable. Wind is often strongest in one or two months out of the year. The 
current index averages over all months. This is potentially a source of error, but 
because the index used is consistently above the threshold for sand transport it is 
assumed to be a small source of error. As long as wind speed is sufficiently high 
the amount of vegetation cover will dominate sand transport. 
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Figure 5.2: Comparison of Mean Monthly Maximum wind speed index to the top 
2% of wind speeds 


5.5.2 Sparse Spatial Distribution 

This study required spatial wind maps; the only available data are point 
specific, thus we needed to spatially interpolate the existing data to create a map 
of the entire wind field. For each scene any wind values within 100km of the image 
boundary were collected and krigged to a 300m grid with IDL’s Krig2D function. 
Despite the low resolution of the grid there was very little change from cell to cell, 
so 300m was deemed to be a high enough resolution for this study. In addition 
each point used appears to be highly correlated with its nearest neighbors thus 
the limited number of available points was not deemed to be a large source of 


error. 
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5.5.3 Homogeneity of Windfield 

The number of data points available to generate a windmap are probably 
adequate because the resulting map is fairly homogenous (fig 5.3). If the data 
points available changed dramatically fi-om one location to the next there would 
be reason to suspect that all of the variability was not captured in the current 
data collection. While it is still possible that this is the case, it appears likely 
that the wind field does indeed vary smoothly from wind station to wind station. 
In addition, the total variability across the High Plains appears to be low relative 
to the variability in vegetation cover. The mean monthly maximum wind speed 
in almost all cases is much higher than the threshold velocity for sand transport 
on bare soil, thus vegetation cover probably dominates the spatial variability and 
any lost variability in the wind field is probably small enough to be negligible for 
this study. 

5.6 Scmd Transport Equation 

One disadvantage to this model is its reliance on a fieldplot-scale sand tran- 
port equation. This equation has limited accuracy and may not represent large 
scale reactivation properly. Effects such as differences in upwind fetch, sources of 
sediment, and downwind sinks which halt the migration of sand are not accounted 
for. Larger scale estimations, such as historical perspectives, implicitly account 
for such effects. 

The sand transport equation used in this study comes from localized mea- 
surements of sand transport and may not be indicative of large scale sand dune 
reactivation. The sand transport equation comes from a study based on individ- 
ual measurements of sand transport with respect to wind and vegetation cover 
on flat ground. Because these measurements do not explicitly account for larger 
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Figure 5.3: The Mean Monthly Maximum wind field across the High Plains ex- 
hibits only minor variability. Here black represents a wind speed of 13 m/s and 
dark red represents 17m/s as measured 10 meters above the ground, a total range 
of 4m/s. Contour lines represent approximately 0.15m/s increments. 
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scale terrain features it is unclear how well they will scale up to estimates of dune 
field activation. Upwind fetch and sediment sources/sinks are important variables 
in large scale reactivation not accounted for by small scale estimates. Bagnold 
(1941) and Shao et al. (1993) showed that sediment erodes much more readily 
when bombarded by sediment entrained upwind, especially fine grained sediment. 
In addition Van Boxel et al. (1999) showed that wind speed varies substantially 
at the surface of a dune due to speed up over the dune crest. Thus dune crests 
may actually activate at a slower regional wind speed due to higher localized wind 
speeds. Using a DEM it may be possible to include for this effect. Other studies 
have investigated larger scale dune reactivation, but cannot account for Landsat 
scale variability and thus are not useful for this study. 

5.6.1 Limited Soil Information 

This study has used only minimal soil parameterization. It may be possible 
to incorporate more detail from the STATSGO soils database. The sand trans- 
port equation developed by Bagnold (1941) allows for calculation of the threshold 
velocity on a bare sandy medium based on mean grain size. Beyond a simple 
mean grain size, the grain size distribution can be important in determining sand 
transport levels. Gillette and Chen (1999) showed that grain size distribution 
within a soil has a major impact on entrainment by wind. This is especially the 
case for dust transport, where the threshold velocity increases with decreasing 
grainsize. In this case larger grains will begin saltating at lower wind speeds and 
their impact can cause erosion of smaller grains which would have not have been 
eroded if the soil was only made up of small grains. At the moment the soil data 
available does not permit the inclusion of this detail, nor are there well developed 
sand transport equations that are a function of grain size distribution. 
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5.6.2 Proper Vegetative Measurements 

It is impossible to make a single measurement of vegetation that describes 
all of the potential effects of vegetation on sand transport. Landsat measurements 
are used to estimate the percent of each pixel which is covered with vegetation, 
but this is not necessarily what is important in the estimation of sand transport. 
If within a single pixel twenty percent of the area is covered with a low standing, 
but dense, grass, large portions of that pixel may still be vulnerable to eolian 
erosion. However, if twenty percent of another pixel is covered with taller grasses 
and shrubs, those plants will have a substantial effect on windspeeds nearby and 
thus the remaining eighty percent of the pixel, though not covered by vegetation, 
will be less likely to experience substantial eolian erosion. Thus many studies of 
sand transport use a measure of vegetation described as the lateral cover. This 
measure incorporates the average height of vegetation as well as the areal coverage 
to provide a better estimate of the amount of vegetation that is actually blocking 
sand transport. 

While it is impossible to measure vegetation height directly using Landsat 
data, it may be possible to estimate this by several other means. Different lo- 
cations through out the High Plains can be characterized by different vegetation 
types. It may be possible to associate an efficiency factor with each area and 
multiply this by the percent cover as measured by Landsat to develop a better 
estimate of lateral cover. In addition a better spatial estimate of this efficiency 
factor may be possible by utilizing multi-angle data from the MISR sensor to es- 
timate structural parameters, or by measuring the phenology of vegetation cover 
at a moderate resolution with the MODIS sensor to determine vegetation type at 
a moderate spatial resolution. 



Chapter 6 


Conclusions 


6.1 Physics based modeling potential 

This goal of this study was to develop and analyze a model for potential 
sand transport in the High Plains of North America. This study has illustrated 
an attempt to bring together the necessary components to run a simple, physics 
based model to estimate future sand transport across the High Plains. This study 
has shown that better modeling of future vegetation is a subject that needs sub- 
stantial work, and has suggested several methods by which this modeling might be 
improved ; a more detailed incorporation of the spatial variability in vegetation 
and soil type, higher temporal resolution satellite measurements and model esti- 
mates, and incorporation of recent rainfall data when estimating vegetation cover 
from NDVI. In addition, further improvements could be made to the model as a 
whole with a better understanding of the effects of vegetation on sand transport 
and by integrating sand transport over several years of wind data. 

6.2 Spatial resolution 

Landsat scale spatial resolution is required to study sand dunes in the High 
Plains. Small, localized dunefields such as those in Eastern Colorado effectively 
disappear with larger (>500m) pixel sizes. Dune crests may be the only areas 
that stand out. If these areas were averaged together with surrounding, vegetated 
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areas by low spatial resolution satellite imagery, the true instability of these dunes 
would not be captured. In addition it is not possible to map land use with enough 
accuracy using lower spatial resolution imagery because even large pivot irrigation 
systems are too small to detect. However, it may be possible to incorporate high 
temporal, low spatial resolution with the high spatial resolution data from Landsat 
to further improve the model. 

6.3 Temporal Resolution 

The benefit of increasing satellite temporal resolution is unclear. The simple 
test performed here showed no relationship between the number of images used 
and the correlation coefficient when regressing Century and NDVI values. It is 
possible that extremely high frequency imagery, such as bi-weekly MODIS data 
would be a substantial benefit to the model, but it is clear that simply collecting 
twice as many Landsat images does not substantially improve the results of the 
present model. 

6.4 Areas of instability 

This model has also mapped areas of potentially higher instability, though 
the present map is not suitable for planning purposes. At the moment careful 
interpretation of the output is required because of the poor correlation between 
NDVI and Century values. Some regions of apparent instability may only appear 
that way due to recent weather conditions or local soil color. It appears that 
the largest areas at risk are in the Nebraska Sand Hills, and portions of Eastern 
Colorado and Western Texas. Further areas at high risk are spread out across 
the High Plains and generally highly localized. This model does not map areas 
that are currently farmed, thus large areas may be at risk that are not shown in 
this model. In addition this model only maps areas of potential sand transport. 
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Dust, which has important implications for human health and well being, is not 
mapped by the present model. 

6.5 Future Work 

One of the major results of this thesis is the understanding of our current 
limitations in this field, and the suggestion of possible methods of improvement. 
There are many holes in our current understanding and in the available data that 
need to be filled before an accurate assessment of dune stability can be made. 
The area that needs the most work is the prediction of future vegetation cover. 
In addition, a better method of incorporating a time series of wind measurements 
and a denser network of weather stations would be beneficial. 

The prediction of future vegetation could probably be improved with a more 
accurate model of vegetation response to climate, and a higher temporal resolu- 
tion sequence of satellite and model estimates. One of the first improvements 
that needs to be investigated is the inclusion of high temporal resolution MODIS 
data. Though this would mean decreasing the spatial resolution of the model, but 
the increased ability to estimate the relationship between climate and vegetation 
would substantially out weigh this loss. 

A better estimate of wind speed over time should be used. The current 
model only includes an index of windspeed. Winds may be highly seasonal in 
nature, so an estimate for each month should be used. Work needs to be done 
to determine the best wind index for this type of study, and what the effects of 
seasonal variation are. 
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